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Abstract. Radiative shock waves show a strong cooling 
instability at temperatures above approximately 2 • lO^K. 
We numerically investigate this instability by simulating 
different astronomical objects in which colliding flows play 
an outstanding role: Wind bubbles, supernova remnants, 
and colliding winds. Computing the flow of each object 
over a large part of its evolutionary time and resolving all 
physically relevant scales, we find several phenomenologi- 
cally different types of this instability. If two smooth flows 
collide, the instability follows a periodic limit cycle with 
two modes being important. The connection between the 
radiative loss function and the mode and type of the over- 
stability is discussed. The collision of non-smooth flows 
can temporarily result in an aperiodic evolution of the 
system. After a characteristic relaxation time the insta- 
bility then becomes periodic again. Such disturbances as 
well as violent types of the instability can excite oscilla- 
tions of the thin layer of cold compressed gas downstream 
of the shock, which in turn can influence the stability of 
the radiative shock. 

Key words: Hydrodynamics - Instabilities - Shock waves 
- ISM:kinematics and dynamics - Numerical methods 



1. Introduction 

1.1. Radiative shocks 

Shock waves become radiative if the dynamical time of 
the flow exceeds the cooling time of the shocked matter. 
Radiative shocks are present in a variety of astrophysi- 
cal objects: e.g. in supernova remnants, in colliding winds 
of binary star systems, in accretion processes, in all wind 
driven structures, in young stellar objects, in galactic and 
extragalactic jets, in star burst galaxies, in active galac- 
tic nuclei. Usually, such shocks are sufficiently strong to 
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substantially ionize and heat the matter. As the heated 
matter cools, it is compressed into thin dense layers. The 
radiative cooling itself and the presence of a dense layer 
of cooled gas greatly influence the emitted spectrum. Ra- 
diative cooling of shock heated matter leads to emission 
mostly in the frequency range between X-ray and UV. 
The cold layer emits mostly in the radio, IR, optical, and 
UV. Therefore, a better understanding of radiative shock 
waves, therefore, is crucial for the interpretation of the 
observed spectra. 

Of similar significance is the influence of radiative 
shock waves on the dynamical shaping of their environ- 
ment. The cold layer is much denser than the surrounding 
matter and often prey to instabilities of various types. In 
addition, the neighboring shock heated zones are often 
subject to thermal instabilities. Such effects can be cru- 
cial for the distribution and the state of the matter in the 
vicinity of radiative shocks. 

In this paper, we study the thermal instability of ra- 
diative shock waves in the context of colliding flows. Ra- 
diative cooling is included as e = 7V^A(T), where A(T) is 
piecewise linear on a log-log scale. The necessary spatial 
resolution to resolve all the relevant length and time scales 
is provided by an adaptive mesh refinement algorithm. In 
contrast to previous investigations, we follow the long term 
evolution of the entire interaction zone. In the case of a 
decelerating interaction zone this leads to an automatic 
scanning of wide ranges of the steady state shock velocity. 
This enables us to make a systematic study of the cool- 
ing instability. It also allows us to investigate the mutual 
influence of the cooling layer and the cold dense layer. In 
agreement with other authors, we find that for smooth 
flows the system is overstable: the instability is oscillatory 
and two different modes are important. The phenomenol- 
ogy of the overstability, however, may be very different 
within one mode which leads us to the introduction of phe- 
nomenological types. In the case of density disturbances 
upstream of the interaction zone, a temporarily aperiodic 
evolution can result for which we give constraints. Finally, 
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we emphasize the importance of the cold dense layer for 
the thermal instability of the cooling layer. 



1.2. The thermal cooling instability 

Catastrophic cooling and the ovcrstability of radiative 
shock waves due to thermal instability of radiative cool- 
ing have been known for several decades. First numerical 
investigations of catastrophic cooling in shock heated plas- 
mas were carried out by Falle (1975,1981). First numerical 
investigations of the overstability of radiative shock waves 
were made by Langer et al. (1981,1982) for the case of 
accretion flows onto compact objects. 

Chevalier & Imamura (1982) (CI in the following) per- 
formed a linear stability analysis for a planar flow against 
a wall as a model for accretion flows. They applied a cool- 
ing law of the form e = ■ T'^, and found the cool- 
ing layer between the accretion shock and the wall to 
be thermally unstable for (3 < j3c. An infinite number 
of modes with successively higher eigenfrequencies were 
found, among them a fundamental and a first overtone 
mode, for which (3^ ^ 0.4 and 0.8 respectively. Numerical 
simulations carried out by Imamura et al. (1984) were in 
agreement with these results. Wolff et al. (1989) and Ima- 
mura & Wolff (1990) extended these investigations to two- 
temperature models with heat conduction and different 
cooling processes such as Compton cooling and relativis- 
tic Bremsstrahlung. They still found the same modes of 
the overstable oscillation of the radiative accretion shock. 

In a beautiful analytical paper, Bertschinger (1986) 
extended the analysis to interstellar shocks. He showed 
that under certain conditions there exists a piecewise self- 
similar solution in spherical symmetry. The flow pattern is 
more complex than for accretion flows: a leading shock, a 
cooling layer behind the shock, a hot or cold interior, and 
a thin layer of already cooled, compressed gas between the 
interior and the cooling layer. Performing a linear stability 
analysis on this solution that allowed for disturbances in 
all three spatial directions, he found two types of insta- 
bilities, a thermal and a dynamical one. For the thermal 
one he found the same unstable modes in the radial di- 
rection as Chevalier & Imamura (1982). Surprisingly, he 
found that no new modes are introduced by nonradial dis- 
turbances. But modes stable to radial perturbations can 
be unstable to nonradial perturbations. 

On the numerical side, significant progress has been 
made by including detailed radiation processes into the 
analysis. Innes et al. (1987) (IGF in the following) and 
Gaetz ct al. (1988) (GEC in the following) considered 
ionization and recombination processes as well as time- 
dependent non-equilibrium cooling. They agreed on the 
onset of the instability for shocks having velocities above 
130 km/s - 150 km/s. GEC found this instability to be os- 
cillatory which is in agreement with further investigations 
presented by Innes (1988). 



Recently, Plewa (1995) studied strong shock waves 

{vs > 1000 km/s) running into dense material and found 
overstable behavior also under these conditions. 

In spherical symmetry, stability analyses of radiative 
shocks were performed by Imamura & Wolff (1990), Houck 
& Chevalier (1992), and Dgani et al. (1995). Spherically 
symmetric cooling shells have significantly different sta- 
bility properties as long as their size is comparable to the 
relevant radius. For small cooling shells,the stability be- 
havior approaches the planar limit behavior. Attempts to 
analyze the thermal cooling instability of two-dimensional 
planar shocks have been made by Strickland & Blondin 
(1995). 

The organization of the article is as follows: In Sect. 2 
we briefly describe our physical model. In Sect. 3, this is 
followed by some notes on the applied numerical methods. 
In Sect. 4 the basic flow patterns are described and the 
sample of computations is presented. The thermal cool- 
ing instability for the case of smooth flows is analyzed in 
Sect. 5, and in Sect. 6 the same is done for the case of 
disturbed flows. The evolution of the cold dense layer is 
investigated in Sect. 7. In Sect. 8 we discuss our results, 
and, finally, we present the conclusions in Sect. 9. 

2. The physical model 

2.1. The governing equations 

Our physical model is based on the Euler equations and 
the use of parameterized radiative loss functions (RLFs in 
the following). We numerically solve the Euler equations 
in spherical symmetry: 



dp ^ dpVr 
dt dr 
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-e{T,N) + q{T). 



(1) 
(2) 
(3) 



Here, p denotes the mass density, w,- the radial velocity, 
p the pressure, e{T,N) the cooling term, and q{T) the 
heating term. The last two are discussed below. The total 
energy density E is the sum of the thermal energy density 
U and the kinetic energy density: 

E = U + pv^/2. (4) 

p and the number density N are connected by 

p = NuruK, (5) 

where mn denotes the mass of the hydrogen atom and p 
the mean atomic weight. In our calculations we always use 
/i= 1.3. 
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This set is closed by a polytropic equation of state 
(with 7 = 5/3 in all our calculations), 



Radiative Loss Function 



u = 



p 



7-1 



The temperature is given by the ideal gas law as 

Pfj,mu 



(6) 



(7) 



This model excludes heat conduction which is probably 
important for strongly shocked astrophysical plasmae. 

2.2. The cooling model 

We describe the radiated energy density per time, e, by 
a temperature dependent RLF A(T) times the number 
density squared. 



(8) 



For the RLF A(T) we use the approximation 

A(T) = Ko,{T)TP\ T g [Ti,Ti+i]. (9) 

RLFs have been published by several authors. The func- 
tions vary significantly depending on the different phys- 
ical processes which are included. The assumed elemen- 
tal abundances are of great importance for the RLF. The 
first radiative loss function (RLFl) we used is a fit to an 
optically thin, time independent RLF taken from Cook 
et al. (1989) which uses solar photospheric abundances. 
The second loss function (RLF2) is a fit to an optically 
thin, time dependent radiative loss function taken from 
Schmutzler & Tscharnuter (1993) which uses Allen's abun- 
dances. This RLF has been obtained by a thermally self- 
consistent, time-dependent computation of the cooling of 
a plasma from 10^ K to 10^ K. Both RLFs are shown in 
Fig. 1 and their exact values are given in Tables 2 and 3 
in the appendix. In all our calculations e is set to zero for 
temperatures below a threshold temperature Tc- 

Our cooling model neglects processes which cool only 
linearly in the density, e.g. cooling by inverse Compton 
scattering. The use of RLFs also completely neglects ion- 
ization and recombination processes. Therefore, we prob- 
ably overestimate the post-shock temperature as the ki- 
netic energy converted by the shock is only used to heat 
the matter. Moreover, in reality, the newly shocked mat- 
ter is likely to be underionized whereas gas on its cooling 
track tends to be overionized. Consequently, the gas may 
radiate at a considerably different rate from the equilib- 
rium value used for RLFl. It probably also radiates at 
a considerably different rate than computed by the time- 
dependent RLF2, since the two cooling histories and thus 
the ionization structures do not correspond. The conse- 
quences of these simplifications are discussed in Sect. 8. 



log Temperoture [Kelvin] 



Fig. 1. The two RLFs used: Cook et al. (1989) (solid line) and 
Schmutzler and Tscharnuter (1993) (dashed line). 



Therefore, our physical model is more approximate 
than the models of IGF and GEC. However, the use of 
an RLF has the big advantage to be simple enough to al- 
low an analysis of the connection between the RLF and 
the response of the system. At the same time, it exhibits 
already most of the features which are characteristic for 
radiative cooling processes behind shocks. Moreover, the 
computational costs are sufficiently low to allow a system- 
atic study. 

For later use we give the expression for the cooling 
time T (by use of eqs. 8, 9 and U = 3/2fcT), 



t{T,N) = 



T 



3fc Ti-/^' 



\dT/dt\ 2Ao, N 



(10) 



2.3. The heating model 



To account for radiative heating of the flow we have in- 
troduced a heating term which, for temperatures below 
a certain threshold temperature Th < Tc, sets the gas 
temperature to Tjj. 

3. The numerical methods 

3.1. The AMRCART code 

Equations 1-4 and 6 are solved by the spherically sym- 
metric version of our AMRCART code^ on the basis of 
a Cartesian discretization. AMRCART combines a mod- 
ern high resolution Eulerian MUSCL scheme developed 
by Colella & Glaz (1984) with the adaptive mesh refine- 
ment algorithm by Berger (1985). It automatically refines 
the spatial grid and the time-step when the local trunca- 
tion error is estimated to be bigger than a given threshold. 
This code was originally designed by Berger & LeVeque 
(1989) for pure gas-dynamics in two space dimensions. It 
was adapted by Walder (1993) to compute complex astro- 
physical problems, including reactive flows and arbitrary 
physical source terms in one, two, and three space dimen- 
sions. 



^ A public version of the code will soon be available. Requests 
please to walder@astro.phys.ethz.ch . 
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Fig. 2. Illustration of the mesh refinement algorithm. Shown 

arc a computed density profile and the numerical mesh. Three 
levels (out of five) of refinement are shown. The mesh size is 
automatically refined by the code. 

The mesh refinement algorithm is crucial to the quality 
of our simulations. Fig. 2 shows how finer meshes are cre- 
ated to improve the spatial resolution where it is needed. 
In regions with refined meshes the time-step is reduced 
by the same factor as well. Thus, each time-scale which 
is dynamically important can be resolved numerically. We 
emphasize that even in regions with extreme high densi- 
ties, and at temperatures at which the RLF A(T) is at 
its maximum, the radiative energy loss of each cell during 
one time step At is at most one percent of the current 
thermal energy in this cell in all our calculations. 

3.2. Some technical details 

In all presented results the center of symmetry is located 
outside the left boimdary of the computational domain. At 
the left boundary we use reflecting (supernova remnant) 
or inflow boundary conditions (wind bubble and collid- 
ing winds). We use free outflow conditions at the right 
boundary. 

3.3. A note on the numerical error 

Numerically obtained solutions are approximate solutions. 

Therefore, one should always convince oneself that the 
numerical solution is sufficiently close to the correct one. 
Figure 3 shows the significant, even qualitative difference 
between a high- and a low-resolution numerical solution of 
an unstable wind bubble (model WBl of Table 4). They 
have been calculated with a finest spatial discretization of 
about 2.4 • 10^^ cm and 1.92 • 10^"' cm respectively. The 
size of the computational domain is 10^" cm. In both cal- 
culations the cfl-number has been fixed to 0.2. Thus, the 
time-steps differ by a factor of 8 as well. Note that even 
the low-resolution calculation has still a remarkably good 
resolution compared to most computations given in the 
literature. 

This large error can be pinned down to the numeri- 
cal smearing of the boundary interface between the cool- 
ing layer and the cold compressed layer. In connection 




yeo rs 



Fig. 3. High- (top) and low-resolution (bottom) simulation 
of the wind bubble model WBl. Shown is the computed 
size of of the radiative shock. The high-resolution simulation 
shows a bigger oscillation amplitude and the oscillation lasts 
nearly twice as long as in the low-resolution simulation. Note 
that in the high-resolution case the oscillation continues after 
IIO'OOO years but is substantially smaller after having under- 
gone a mode- transition (see Fig. 7). 

with stiff source terms (radiative cooling in our case), a 
too broad smearing of the interfaces can cause unphysical 
wave speeds (see e.g. LeVeque & Yee 1990, Colella et al. 
1986 or Klingenstein 1994). This problem is inherent to 
each numerical algorithm which is based on interface cap- 
turing by finite volumes. 

We estimate the errors of our calculated interface 
speeds to be at most a few percent. 

4. Investigated flows 

Simulations of wind bubbles (WB), supernova remnants 
(SNR) and colliding winds (CW) are the basis for our 
study. These flows are well documented in the literature 
and we give here only the characteristics needed for this 
work. A WB (see Fig. 4) is the structure which results 
when a strong stellar wind blows into the uniform inter- 
stellar medium. Similarly, a SNR is the result of a strong 
explosion driving expelled gas against the uniform inter- 
stellar medium. In the case of CW, a fast wind sweeps 
up a slower predecessor wind from the same star. In the 
adiabatic case, all three flows are self-similar. The leading 
shock decelerates in the WB-casc according to Vs oc 
(vs denotes the shock velocity, t the time) and in the SNR- 
case according to Vg oc t~^/^. Ryu & Vishniac (1991) give 
a nice overview of these flows. In the CW-case, the ve- 
locity of the leading shock is constant (Chevalier & Ima- 
mura 1983). If cooling is included the flows are no longer 
self-similar. However, Bertschinger (1986) shows that un- 
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Fig. 4. Leading part of the interaction zone in the radiative 
case. For details see text. 



der certain conditions radiative SNR and WB can be de- 
scribed by piecewise self-similar solutions. 

Wc arc interested in flows where the cooling layer is 
small against the radius of the structure. In this case - de- 
spite their different driving mechanisms - WB, SNR and 
CW have a similar leading part of the interaction zone 
where the two flows collide. The leading part of the inter- 
action zone is illustrated in Fig. 4 at the example of the 
WB. In all three cases, it consists of a leading shock (LS), 
a hot leading shell in which the shocked material cools 
(HLS) and a 'cold' dense layer with compressed, already 
cooled shocked material (CDL). The HLS is subject to the 
cooling instability, the CDL can play an important role for 
the evolution of this instability. We are interested in these 
two parts. The contact discontinuity (C) marks the rear 
end of the leading part of the interaction zone. Notice that 
in our study the reverse shock is nearly adiabatic since - 
due to the much lower density and the much higher tem- 
perature - the cooling time is considerably longer than 
the dynamical time. In the following we do not discuss 
this reverse, driving part of the flow. 

We list the exact settings of all our calculations in Ta- 
ble 2 in the appendix. We performed SNR and WB sim- 
ulations with different parameters up to the time where 
the thermal instability ceased (see Fig. 3). CW were per- 
formed over a dynamical time which is equal to ap- 
proximately 600 cooling times of the leading shell. All 
three flows have been calculated with both RLFs given 
in Sect. 2.2 and in the appendix. The results presented 
below are the quintessence of all these calculations. Of ma- 
jor importance is the WB calculation. Since the structure 
moderately decelerates, it is possible to scan a wide range 
of post-shock tcmperatines (in our case from T « 10^ ''' K 
to T « 10^ K) while the global conditions remain nearly 
constant over one oscillation period of the instability. 

5. The cooling instability of the leading shell for 
colliding smooth flows 

In this section we present the results for the case where two 
smooth flows collide. In Sect. 6 we investigate the insta- 
bility for disturbed flows. At present, smoothness means 



that during one oscillation cycle the conditions for the ra- 
diative shock remain nearly constant. In Sect. 6 we give 
more thorough constraints. 

In agreement with earlier investigations, our system- 
atic study shows that the cooling instability for smooth 
flows - if present - manifests itself as an overstability: All 
quantities, e.g. the size of the cooling layer or the radi- 
ated energy, evolve in a periodic way with fixed period 
and amplitude. The linear stability analyses of CI and 
Bertschinger (1986) show that the period of the oversta- 
bility of the radiative shock is determined by dimension- 
less quantities 6i {i = 0,1,2, , ...) which are eigenvalues of 
the boundary value problem describing the flow between 
the radiative shock and the rear end of its cooling layer. 
One usually speaks of the fundamental mode (F-mode, 
i = 0) and of overtone modes (10-mode, 20-mode, etc). 
The governing eigenvalue 5 is related to the steady state 
position of the shock, xq, the velocity of the inflowing mat- 
ter relative to the shock, Ui„, and the oscillation period, T, 
by the relation S = xo/uin ■2n/T. Previous investigations 
have shown that even in the non-linear regime the period 
is in good agreement with the fundamental (F) and the 
first overtone (10) mode determined by linear stability 
analysis. 

However, our analysis shows - also in agreement with 
earlier investigations that the overstability can show dif- 
ferent phenomenological manifestation when oscillating in 
one particular mode. Therefore, we propose to introduce 
an additional classification of the overstability, namely the 
notation of types. 

5.1. Different types of oscillations 

Following phenomenological criteria, we propose five dif- 
ferent oscillation types. The most striking differences be- 
tween the different types are: the oscillation amplitude, 
which we express as A ={qmax-(lmin) /imin where q is a 
characteristic quantity of the cooling layer, whether the 
oscillations are sinusoidal or not, and the oscillation mech- 
anism. Fig. 5 and Fig. 6 illustrate the characteristics of the 
different types. In the following, we give only a brief de- 
scription of some aspects of these types and refer to earlier 
authors (see Sect. 1.2) for more detailed descriptions. 

The catastrophic C-type is characterized by a dramatic 
collapse of the temperature at a certain moment of the os- 
cillation cycle. In wide ranges of the cooling layer the tem- 
perature drops nearly to nebular temperature on a very 
short time-scale. The density is hardly affected despite the 
rapid change in pressure due to the strong cooling (snap- 
shots 4 and 5). A huge pressure cavity evolves and a sec- 
ondary shock forms at the boundary to the cold layer as 
the pressure equilibrium starts to be re-established (visible 
in snapshots 3 - 6) . The renewed growth sets in when the 
secondary shock in front of the cold layer hits the undis- 
turbed medium again (snapshot 1). The amplitude A is 
large and can reach up to two orders of magnitude. None 
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snapshot 1 snapshot 2 snapshot 3 
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Fig. 5. a) Timo-sorics over one oscillation cycle (from left to right) for each different type of the overstability. The series 
continue on the next page. Shown are the density- (solid line) and temperature-profiles (dashed line) of the cooling layer. From 
top to bottom: Sl-type (F-mode with superimposed lO-mode), I-type (F-mode), C-type (F-mode), M-type (lO-mode), and 
S2-type (lO-mode). The series correspond to the cycles shown in Fig. 6. They start approximately at the minimum extent of 
the shell (phase </> = 0) and end at the next minimum {(j) = 1). (Caption continues next page.) 
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Fig. 5. b) Continuation of Fig. 5. a) of previous page. The length of the cooling layers are measured in units of the lengths 
of their first minima {<f) = 0). Notice the large difference in amplitude between the different types of oscillation for the F-mode. 
For the lO-mode, the amplitude is generally small. Note that the cooling layers may have different sizes at phases ^ = and 
= 1. This is due to the influence of the dynamics in the cold dense layer (see Sect. 7) and the general deceleration of the 
structure. Although the types arc independent of the particular flows we give in brackets the model from which we have taken 
the example: Sl-type (CW), I-type (WBl), C-type (WBl), M-type (WB2f) and S2-type (WBl). 
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Fig. 6. The different types and modes of the cooling overstability shown in some characteristic quantities. The graphs correspond 

to the cycles shown in Fig. 5, From top to bottom: Size of the leading shell, radiated power (of the entire leading shell), post shock 
temperature, shock-velocity, average density of the leading shell and mass-flux from the leading shell into the cold compressed 
gas layer. Notice the highly different amplitudes of the oscillations in the different types. 
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of the quantities oscillates sinusoidally. In particular, the 

radiated energy is rather peaked. For the RLFl we find a 
C-type oscillation for stationary post shock temperatures 
rf/ > 4 • 10^ K. The upper limit probably lies higher than 
Tg^ K, 4 • 10® K, the highest limit post shock temperature 
we have investigated for RLFl. For the RLF2 the C-type 
is present between 3 • 10^ K < Tf/ < 5 • 10^ K. We observe 
the C-type only in the F-mode. 

As another extreme, we observe the smooth S-types. 
There, the shell remains hot during the entire oscillation 
cycle and no pressure cavity evolves. No secondary shocks 
are observed. All quantities oscillate rather sinusoidally. 
The Sl-type can be associated with the F-mode, the S2- 
type with the 10-mode. A is always of order one, but 
is generally smaller for the S2-type. In the case of the 

51- type, as the shell shrinks, a pressure gradient evolves 
between the cold layer and the leading shock. The force 
it exerts leads to the renewed acceleration of the leading 
shock. In the case of the S2-type, on the other hand, a pres- 
siire wave evolves at the boundary to the cold layer which 
later hits the leading shock and accelerates it. This leads 
to the phase-shift between the energy loss and the size of 
the shell (sec Fig. 6). The differences between the Sl-type 
(F-mode) and the S2-type (10-mode) correspond to the 
results of the linear analysis of CI which predict knots and 
phase-shifts for overtone modes. In the RLFl we find the 

52- typefor temperatures 1.7- 10^ K < Tf/ < 2.5-105K. 
Note that the stability limit is determined by the insuffi- 
cient spatial resolution in the simulation and may be some- 
what lower in reality. We find no Sl-type in the RLFl. For 
the RLF2, we find the Sl-type for temperatures between 
1.8 • 10^ K < Tf/ < 2.2 • 10^ K and the S2-type for tem- 
peratures 6.6 • 10^ K < Tf/ < 8.5 • 10^ K. 

The intermediate I-type occurs during the transition 
from C-type to S-types. The leading shell always remains 
rather hot, but - as in the C-type - the shell vanishes apart 
from the secondary shock. The amplitude of the oscillation 
is still rather big. We find the I-type for temperatures in 
the range of 2.5 • 10^ K < Tf/ < 4 • 10^ K for the RLFl 
and in the range of 2.2 • 10^ K < Tf/ < 3 • 10^ K for the 
RLF2. Also the I-type is observed only in connection with 
the F-mode. 

In the mixed M-type the leading shell is divided into 
two parts, namely a high density part near the cold layer 
and a low density part behind the leading shock. While the 
low density part remains hot the high density part strongly 
cools (snapshots 2 and 3). When this cooled dense mat- 
ter is swallowed by the cold dense layer a pressure wave 
or a weak shock is generated (snapshot 4) which travels 
through the leading shell towards the leading shock (snap- 
shots 4, 5 and 6). On its way it enhances the density of the 
areas it passes. Finally, it hits the leading shock and accel- 
erates it (snapshot 6). By now the shell is again divided 
into a high density and a low density part (snapshot 6 
and 1). The oscillation amplitude is of order one. Most of 
the energy is released during a very short time. We find 



the M-type for the temperature range between 5 • 10^ K 
< Tf/ < 6.6 • 10^ K of the RLF2, we do not observe it 
for the RLFl. The M-type oscillates in the 10-mode. Like 
the S2-type, the M-type shows knots and phase-shifts. 

The phenomenology of the overstability and, in par- 
ticular, the amplitudes of the oscillations are essentially 
determined by the non-linear terms. Each of the modes 
lives mainly in two different forms, a weak and a strong 
one. The weak form of the F-mode corresponds to the Sl- 
type, the strong form to the C-type, while the weak form 
of the 10-mode to the S2-type and the strong form of the 
10-mode to the M-type. 

5.2. What determines the modes and the types? 

We do not yet have a final answer to this question. How- 
ever, we present some observations. They are mostly based 
on mode-transitions which occur in WB-simulations when 
the post-shock temperature slowly decreases due to the de- 
celeration of the entire interaction region. We begin with 
two general observations. 

• 1) According to our experience, for a given radiative 
loss function, the type of the overstability is fixed by the 
post-shock temperature Tf/ of the (unstable) stationary 
solution alone: Different forms of excitation of the insta- 
bility at the same Tf/ result - after a certain relaxation 
time - in the same mode and the same type. We detect no 
sign of multiple solutions or a, bifurcaMon of the solution. 

• 2) Linear theory predicts that whenever the F-mode 
is present the 10-mode is present as well. However, in the 
strong forms of the F-mode we do not observe any sign 
of a superimposed 10-mode. Only at the transition from 
the I-type to the S-type a slight trace of a superimposed 
10-mode shows up. At this point, we observe no phase 
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Fig. 7. Mode-transition in WBl from the F-mode (I-type) 
to the lO-mode (S2-type) for the RLFl. Shown are the 
timc-cvolution of the size of the leading shell (top) and the 
post-shock temperature (bottom). 
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Fig. 8. WB2f: Stability behavior when applying RLF2 for 
1-10® > T^^ > 5•10^ Shown is the size of the leading shell (top) 
and the post-shock temperature (bottom). Above T^^ « 7-10^ 

the (artificially excited) F-mode (C-type) is strongly damped, 
below this temperature a lO-mode is growing. The further evo- 
lution is shown in the next flgure. 

shift between the modes due to the non-rational ratio of 
the eigenfrequencies. The phase shift only sots in when the 
overstability is already very smooth (Fig. 7) . The physical 
reason for the absence of overtone modes in the strong 
forms may be the total destruction of the leading shell 
and, therefore, the history after each oscillation cycle. 

The next two points are dealing with the transition of 
types and modes when scanning RLFl in the temperature 
range where the logarithmic slope changes from a negative 
to a positive value (see Fig. 7). 

• 3) We observe a mode and type-transition from 
the F-mode (I-type) to the 10-mode (S2-type) at « 
2.5 • 10^ K when applying RLFl. That is approximately 
the temperature where (5 jumps from a strongly negative 
value to a moderately positive value. For this transition, 
the change of j3 associated with the post-shock temperature 
seems to he critical. 

• 4) The absence of the F-modc below T^f « 2.5 -10^ K 
for RLFl is surprising. The post-shock temperature varies 
over the entire oscillation cycle between 1.9 ■ 10^ K and 
2.3 • 10'"' K and even most parts of the shell have always 
temperatures above 7.9 • 10'* K. For these temperatures, /? 
is equal to 0.15 while below 7.9-10^ K it becomes strongly 
positive. Linear analysis predicts for /3 = 0.15 the F- and 
the 10-mode to be present. Indeed, when using a RLF 
assuming a constant /3 = 0.15 we observe both modes. We 
conclude that the mode is not governed by the (3 which 
is associated with the stationary post-shock temperature 
alone. The entire cooling layer and, in particular, the rear 
boundary interface may be of similar importance. The 
strongly positive /? for T < 7.9 • 10** K reduces the cooling 
efficiency as compared to /? = 0.15 and in this way seems 
to damp the F-mode. 



Fig. 9. WB2f: Mode-transition from 10-modc (M-type) to 
F-mode (C-type) at Tf/ < 5 - lO'' when applying RLF2. Shown 
is the size of the leading shell (top) and the post-shock tem- 
perature (bottom). (The evolution at earlier times is shown in 
the previous figure.) 

The next three points summarize the transition of 
modes and types when scanning RLF2 in the tempera- 
ture range 1 • 10*^ ^ T^| ^ 4 • 10^ K (see Figs. 8 and 9). 
This is exactly the range where RLF2 forms a 'valley'. For 
T < 4 • 10^ K the logarithmic slope of the RLF2 is equal 
to -1.92, for 4 • 10^ < T < 5.37 • 10^ K it is zero and above 
5.37-10^ K it is equal to 1.12. The planar models, CPF_55 
to CPF_25, prove that the different modes and types ob- 
served in this 'valley region' of the RLF2 are not merely 
transition phenomena but persist if the structure does not 
slow down. 

• 5) The excited F-mode (C-type) at ^ 9 • 10^ K is 

very strongly damped when RLF2 is applied. This suggest 
a stable situation for this temperature, and indeed the 
planar model CPF_55 with T^j = 9 • 10^ K is stable. At 

this temperature we have (3 = 1.12 which is in the stable 
range according to linear theory. However, large parts of 
the cooling layer have temperatures associated with an 
unstable /?. Despite this fact, in this case the (3 associated 
with the post-shock temperature seems to determine the 
stability. 

• 6) Below Til ~ 8.5 • 10^ K a 10-modc of the S2-type 
grows for RLF2 (CPF_45 for T^^ « 8.25 • 10^ K shows 
S2-type behavior). The post shock temperature does not 
leave the region with /? = 1.12 over the entire oscillation 
cycle. Interestingly, nearly in the entire cooling layer the 
temperature is in the range where (3 = 1.12! With de- 
creasing T^j*, the type of the overstability slowly changes 
to the strong form of the 10-mode (M-type) (CPF_35 for 
T^t « 6.1 • 10^ K is clearly of the M-type). (3 associated 
with T^f" is still 1.12! This again indicates that the RLF 
for temperatures below plays a significant role. 

• 7) As can be seen from Fig. 9, at 5 ■ 10^ K 
another mode-transition takes place for RLF2, this time 
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from the 10- (M-type) to the F-mode (C-type). This tem- 
perature corresponds approximately to a change in /? from 
the stable value 1.12 to zero, a value at which - accord- 
ing to linear theory - both modes should be present. The 
change of (3 corresponding to the post-shock temperature 
seems to he important for the mode-transition. 

The emerging picture is not homogeneous. Sometimes, 
the (3 associated with T^f seems to play an outstanding 
role. More often, however, the presence, the mode, and the 
type of the overstability cannot be determined on the basis 
of the local characteristics of the RLF at this temperature. 
We come back to these points in Sect. 8 

6. The cooling instability of the leading shell for 

disturbed flows 

6.1. The set of disturbed flows 

We now discuss the cooling instability for disturbed flows. 
We disturb the flow upstream of the interaction zone be- 
tween position ro and ro -|- Aq by setting the density to 



post shock temperoture 



Pir) = po 



1 + ^ sin TT 



Ao 



if ro <r<ro + Ao (11) 



where po is the density of the undisturbed flow. The distur- 
bance is applied after the system has evolved into a limit 
cycle. The leading shock can hit it at an arbitrary phase 
< (/lo 5J 1 of the oscillation cycle (where = 0,1 arc the 
phases of minimum extent of the leading shell). Note the 
difference to the setting of IGF who let the shock in its 
stationary position run against a, however similar, density 
disturbance. 

We performed a systematic investigation of disturbed 
I- type systems on the example of WBl. All disturbances 
were applied around 69'000 years. The exact settings of 
these models can be taken from Table 5 of the appendix. 
Except for two models (WBl_dO and WBl_dl) the am- 
plitude of the disturbance, ^, is never bigger than 1. The 
ratio of the mass of the disturbance, Mj, to the mean mass 
(over a period) of the hot leading shell, M^g, varies from 
0.017 to 73, the ratio of to the mass of the subsequent 
cold dense layer, M^jj, from 1.4-10"'* to 0.6155. The ratio 
of the length of the disturbance, Ag, to the characteristic 
length of the limit cycle, i — T-v^^, varies from 0.05 to 8.8 
(T denotes the oscillation period and wf^ the stationary 
shock velocity). The shock hits the disturbance at phases 
between (/>o = 0.07 and 0.52. The results presented in this 
section are a summary of the analysis of this set. 

6.2. Noisy, original periodic, aperiodic, and modulated pe- 
riodic evolution 

We found five clearly different responses of the I-type sys- 
tem to disturbances. The time-evolutions of three of them 
are shown in Fig. 10 and a schematic representation of 
all possible responses is given in Fig. 11. If the distur- 
bance is small (in which sense we see below), one ob- 
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Fig. 10. Noisy original periodic, aperiodic and modulated pe- 
riodic evolution of the radiative shock after being disturbed at 
time 68'893 years. Shown are post-shock temperatures of mod- 
els WBl_d4.0.1, WBl_d4.1 and WBl_d6 (from top to bottom). 



serves a rapid, aperiodic, much smaller scaled variation 
of only one cycle of the original oscillation. We call this a 
noisy, original periodic evolution (np). If the disturbance 
is appropriately scaled, the periodicity of the limit cycle 
is completely destroyed and we have an aperiodic evolu- 
tion of the radiative shock (ap). However, after a certain 
time, the system relaxes to the original, periodic evolu- 
tion. Disturbances on a larger length scale and not too 
big amplitude ^ result in a modulation of the limit cycle, 
governed by the slowly changing pre-shock conditions. In 
this case, we speak of a modulated periodic evolution (mp). 
On the other hand, if we have a very massive disturbance 
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Fig. 11. Schematic representation of the different responses 
versus mass and length of the density disturbance. 
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Fig. 12. The aperiodic evolution of the radiative shock and its cooling layer from example of WBl_d4.1. Shown are density- 
(solid line) and temperature-proflles (dashed line) of three successive times. The evolution is governed by the fast partial cooling 
and the waves which are created during this process. Details can be taken from Sect. 6.3. 



on a relatively short length-scale, the leading shell col- 
lapses nearly immediately (co). The history of the leading 
shell is destroyed and the shell behaves according to the 
new global conditions. For example, in model WB_dO the 
accumulation of the mass of the disturbance by the cold 
dense layer leads to a slow down of the entire interaction 
zone and the leading shock becomes stable. Finally, if we 
encounter a small mass disturbance on a large length scale, 
the; systcini shows nearly no reaction (nr). We note already 
here, that big disturbances may significantly disturb the 
cold dense layer downstream of the radiative shock wave 
(see Sect. 7). 

6.3. The aperiodic evolution 

The oscillation mechanisms of the np- and the mp-cases 
are similar to the mechanisms described in Sect. 5. Here 
we describe the aperiodic response of the radiative shock 
at the example of model WBl_d4.1. It evolves very similar 
to what was described by IGF. We concentrate on some 
aspects important to the understanding of the constraints 
for the evolution of disturbed shocks. 

The size of the initial density disturbance is crucial 
for the time evolution illustrated in Figs. 10 and 12. As 
the density disturbance passes through the leading shock 
it shrinks in size while being compressed. Cooling of the 
hot shell starts in the region with enhanced density. In 
fact, cooling is so fast that the flow cannot maintain pres- 
sure balance and a pressure cavity forms in the cooled re- 
gion. A sheet of cold and initially over-compressed matter 
is then created within the hot leading shell (snapshot 1 
of Fig. 12). The relaxation of the over-compressed sheet 
generates two pressure waves or weak shocks traveling 
forwards and backwards in the hot shell. Their strengths 
are a direct measure of the violence of formation of this 
second dense sheet. The dynamical consequences of these 
waves are twofold. Firstly, temperature and density in- 
crease again in those regions where the waves pass. The 
consequence for the cooling scale cannot be estimated eas- 
ily, since in this region of the radiative loss function the in- 
crease in temperature and density are working in opposite 
directions. For RLFl and in this temperature range the 



cooling scale likely is to decrease behind weak waves, it is 
likely to increase behind strong waves. In our example, the 
first possibility applies. A second important consequence 
of these waves is encountered when they eventually hit the 
leading shock. The pushing wave transfers part of its en- 
ergy to the leading shock which is accelerated. As a conse- 
quence of the shock jump conditions, the post-shock tem- 
perature is higher and the density is lower in the freshly 
shocked part than in the older parts of the shell (compare 
snapshot 2). The further development can hardly be fol- 
lowed in detail. The older, denser parts cool faster than the 
newer parts. Fiu'ther cold, compressed sheets are formed, 
usually near the boundary to the cold layer. Their forma- 
tion creates again new waves and the process continues 
(snapshot 3). 

Our investigation proves - for the I-type range and for 
parameterized cooling - that this process does not last 
forever and that the limit cycle finally recovers. 

6.4- The constraints for the evolution and relaxation times 

As can be seen from Table 5 in the appendix, models hav- 
ing a length of the disturbance between 5 and 20 percent 
of the characteristic length of the limit cycle and hav- 
ing masses between 10 and 40 percent of the mass of the 
leading shell react with an aperiodic evolution to the dis- 
turbance (case ap). 

To achieve the most efficient aperiodic evolution, the 
region of enhanced density has to be moderately smaller 
than the hot shell, since otherwise the entire leading shell 
would be affected nearly similarly and we would encounter 
a modulation of the original oscillation (case mp) or a 
collapse of the shell (case co). On the other hand, the 
region with enhanced density has to be large enough to 
be able to create secondary waves as powerful as possible, 
since these waves are driving the aperiodic evolution. 

For the aperiodic evolution, we define the relaxation 
time, Trei, as the difference between the time when the 
aperiodic evolution first sets in and the time when the reg- 
ular periodic evolution starts again, r^ei can be taken from 
Table 5. The biggest ratio of Trei/To (tq being the period 
of the limit cycle) we measured is 3 (model WBl_d4.1). In 
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this case, the ratio of the characteristic lengths, Ao/^c, is 
0.2, the mass ratio, Mj/Mj-^j^^,, is 0.4. Comparing models 
WBl_d3, WBl_d2.2 and WBl.d4.1 we notice the influence 
of the different masses. A larger mass of the disturbance 
leads to the generation of stronger pressure waves, and the 
relaxation time increases. 

7. The cold dense layer (CDL) and the hot leading 
shell (HLS): An interacting system 

Downstream of the radiative shock and its cooling layer, 
a thin layer of cold, compressed gas is established (see 
Fig. 4). In this section we study the dynamics of this layer, 
in particular its interaction with the thermally unstable 
radiative shock. 

7.1. Mass feeding of the CDL 

We define the boundary between the cold dense layer and 
the hot leading shell to be at the surface where cooling 
is balanced by heating. Through this boundary interface, 
there is a net mass-flux from the cooling layer of the ra- 
diative shock into the cold dense layer which is shown in 
the last row of Fig. 6. 

It is interesting that for smooth flows the five different 
types of the overstability manifest themselves also in these 
mass-fluxes (Fig. 6). Most remarkable, there is a very sig- 
nificant difference between the strong and the weak form 
of the overstability. For S-type instabilities the mass-fluxes 
are smooth and vary only moderately over one cycle. In 
the strong forms of the overstability, most of the mass en- 
ters the CDL during a very short time-interval of the cy- 
cle. This is a direct consequence of the sudden collapse of 
large amounts of the hot shell. Again we note a significant 
difference between the modes: In the 10-mode (M-type), 
the peak of the mass-flux is correlated with the largest, in 
the F-mode (I- and C-types) with the smallest extension 
of the cooling layer. Note also that in the I- and C-type 
overstability the onset of the mass-flux at approximately 
one third of the period is accompanied by the formation 
of the secondary shock which is attached to the CDL. Ob- 
viously, the irregularity of the mass flux is enhanced in 
disturbed flows. 

7.2. Response of the CDL to the cooling overstability of 
the leading shell 

The non-steady mass-flux through the interface between 
the HLS and the CDL introduces disturbances into the 
CDL, even when smooth flows collide. Assuming that the 
size of the CDL is much bigger than the HLS, Bertschinger 
(1986) has shown by linear stability analysis that the dis- 
turbances damp on a scale which is much smaller than 
the size of the CDL. When we encounter the same con- 
ditions, our numerical study shows indeed this behavior. 
However, there are many cases where the size of the CDL 
and the HLS are comparable, e.g. in the period after the 
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Fig. 13. Density distribution of the CDL of WBl (I-type in- 
stability) immediately after the leading shell has collapsed and 

the ovcrcomprcsscd shoot of freshly coolod matter is swallowed 
by the CDL (left) and half a period later in the limit cycle 
(right). 



transition of the shock from the adiabatic to the radiative 
regime. Moreover, the analysis in Sect. 7.1 shows that in 
strong forms of the overstability the mass feeding into the 
CDL - and, therefore, the introduced disturbance - is ex- 
tremely pulsed with a large pulse height. When disturbed 
flows collide, these disturbances can become even larger. 
Such disturbances arc clearly beyond the scope of linear 
analysis. Therefore, it is not surprising that under these 
different circumstances the disturbances can survive for a 
significant time. 

7.2.1. Traveling waves and density variation 

As described in Sect. 6.3 the sudden cooling of large parts 
of the hot leading shell leads to the formation of an ini- 
tially overcompressed sheet of cold matter. Apparently, 
the same mechanism also works in the strong forms of the 
limit cycle of the overstability when smooth fiows collide, 
with the difference that the newly cooled, overcompressed 
sheet is in contact with the already present CDL. The 
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Fig. 14. Oscillating size (top) and density-variance a (bot- 
tom) of the CDL. Shown are 18 cycles after the transition 
from the adiabatic to the radiative regime of model WBl. The 
overstability of the leading shell is of the C- and I-type. 
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swallowing of such an overcompressed sheet by the CDL 

is shown on the left of Fig. 13. The thin sheet at the very 
right boundary of the CDL, having approximately twice 
the density of the rest, contains the freshly cooled matter. 
At this time, the density variance in the shell is greatest 
(see Fig. 14). In the following, relatively quiescent phase 
of the limit cycle, the overcompression relaxes by a wave 
traveling through the CDL (right graph of Fig. 13). Even- 
tually, the propagating wave is reflected at the rear in- 
terface of the CDL and interacts with the other waves. 
As a result, the density structure of the CDL is rather 
complicated in the presence of strong forms of the over- 
stability of the radiative shock. The evolution of the den- 
sity variance shown in Fig. 14 exactly mirrors the cycle of 
the overstability of the leading shell. It shows peaks every 
time a new, overcompressed sheet is incorporated into the 
CDL. The figure demonstrates that the density variance is 
rather high as long as we encounter I-type instability. Not 
visible in the figure, the CDL rapidly smoothes after the 
instability changes to S-type. Note also, that we do not 
find a significant difference between the density variance 
of an isothermal CDL (model WBlc) and an adiabatic 
CDL (model WBl). 



7.2.2. Oscillations of the CDL 

The CDL is expected to become compressed when it 
pushes the shock after the collapse. This should lead to 
a considerable shrinking of the size of the CDL, at least 
as long as there is not too much mass in the CDL. In- 
deed, we observe such oscillations (see Fig. 14). However, 
the traveling, reflecting, and interacting waves disturb the 
basic oscillation cycle remarkably. Periods of nearly no 
oscillation are followed by periods of nearly regular oscil- 
lation and vice-versa. When smooth flows collide, we find 
changes in the size of the CDL up to 15 percent as long as 
we encounter I-type instability in the WB and SNR simu- 
lations. Much larger oscillation of the CDL can be excited 
in disturbed fiows as will be shown in the next paragraph. 



1.3. The back- coupling to the cooling layer 

We now investigate how the dynamics of the CDL affect 
the radiative cooling overstability of the hot leading shell. 
Of course, the oscillation of the CDL changes the rear 
boundary conditions of the cooling layer with time, which 
influences the evolution of the radiative shock. This back- 
coupling can be strong in disturbed flows as will be demon- 
strated below. When smooth flows collide, the effect is 
present but never large (See Fig. 3). Only immediately af- 
ter the transition from the adiabatic to the radiative phase 
it is considerable. In this phase the amplitude of the over- 
stable oscillation can change up to 15 percent from one 
period to the next. Although weaker, the effect can be 
observed as long as the overstability is of the I-type. 



7.3.1. The evolution of model WBl_d7 

WBl_d7 is a disturbed model with a modulated periodic 
behavior (see Table 5). The disturbing mass is 7.3 times 
the average mass of the hot leading shell and 6.18 percent 
of the mass of the CDL. The encounter of the interac- 
tion zone with this density disturbance excites an oscilla- 
tion of the CDL with an amplitude of about 15 percent. 
The period of the oscillation approximately corresponds to 
the length of the density disturbance. As demonstrated in 
Fig. 15, this oscillation of the CDL signiflcantly influences 
the further evolution of the cooling overstability. 

Note also the direct consequence of the density distur- 
bance for the overstability (see again Fig. 15). The leading 
shock slows down. The oscillation period becomes smaller 
due to the decreasing cooling time. However, the radiative 
energy release increases due to the density enhancement 
which overcomes the decrease of the emitting volume due 
to the smaller size of the shell. When the interaction zone 
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Fig. 15. Time evolution of the model WBl_d7 after running 

into a disturbance having 6.18 percent of the mass of the CDL. 
1. line: Size of the CDL, 2. line: Velocity of interface between 
CDL and HLS, 3. line: Size of HLS, 4. line: Radiated energy 
per time from the HLS. 
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has passed the density peak, it accelerates again. At 73'500 
years the oscillation of the leading shell has nearly recov- 
ered the state it had before the disturbance. However, the 
absorbed mass has led to a general slow down of the in- 
teraction zone and to the oscillation of the CDL. 

7.3.2. Switch on, switch off system 

The onset of the cooling instability is very sensitive to the 
post-shock temperature. Thus, the instability can switch 
on and off under the influence of an oscillatory cold dense 
layer as we will demonstrate with a CW calculation. At 
the beginning of the evolution shown in Fig. 16 a strong 
wave is traveling backwards and forwards in the CDL. 
Due to the oscillation of the CDL, the velocity of the rear 
boundary interface of the cooling layer is trans-critical: 
If the velocity falls below the critical value, the cooling 
instability stops, if the interface is again pushed over the 
critical value, the instability reappears. The overstability 
is of the Sl-type. The instability immediately becomes 
periodic if the critical post-shock temperature is exceeded. 

The third instable time period is most interesting. The 
instability is excited by a very small overshooting of the 
shell over the stationary value. Although the excitation is 
small, the amplitude of the oscillation is growing up to a 
value which is inherent to the system and the applied cool- 
ing function. Note that the velocity of the rear boundary 
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Fig. 16. The overstability switches on and off under the in- 
fluence of an oscillating CDL. Top: Size of the CDL. Middle: 
Velocity of the interface between the CDL and the HLS. Bot- 
tom: Post-shock temperature of the leading shock. 



of the cooling layer is approximately constant or is even 

decreasing while the amplitude of the oscillation is grow- 
ing. The stationary amplitude of this oscillation is nearly 
the same as the one in the second instable time period. 
The period of the limit cycle increases with time due to 
the decreasing density in the slow wind. 

8. Discussion 

8.1. What determines the types and modes? 

The large sample of computations presented in this pa- 
per strongly indicates: If parameterized loss functions are 
applied to smooth colliding flows, a differently disturbed 
radiative shock wave always relaxes to the same asymp- 
totic solution, whether this solution is stable or overstable. 
However, the time scale of the relaxation may strongly de- 
pend on the size and the sign of the excitation. The sta- 
tionary post-shock temperature T^^ completely fixes the 
stability properties (stability limit, mode, type) of a ra- 
diative shock for a given RLF. 

The cooling overstability in a C-type (I-type) manifests 
itself as runaway cooling. This requires that a significant 
amount of the thermal energy of the cooling layer is ra- 
diated on a time-scale which is short compared to the 
hydrodynamical time-scale. Thus, /? has to be sufficiently 
negative at T^| and in a sufficiently wide temperature 
range below. For example, this condition may be fulfilled 
at temperatures where line cooling of a strong coolant 
just sets in and abruptly enhances the cooling efficiency. 
For these types, the stability properties are completely de- 
termined by the strongly negative logarithmic slope of the 
RLF around T^^. At least, we see hardly any difference 
between RLFl and RLF2 even though the two loss func- 
tions are different for lower temperatures, except for the 
absolute values of the oscillation period which are due to 
the shift in A-direction by half a magnitude. 

Runaway cooling and, therefore, negative /3 are also 
necessary for the occurrence of the M-type. Here, however, 
only the low temperature part at the rear end of the shell 
collapses. To prevent the shell from complete collapse (as 
in the C-type), the cooling time in the hotter front part of 
the shell has to be sufficiently longer than the time needed 
to collapse and reconstruct the rear part. This requires a 
sufficiently positive (3 at temperatures T w Tj*/. These 
conditions may be fulfilled at temperatures between re- 
gions with efficient line cooling for different ions. For this 
type, not only the shape of the RLF at is important for 
the stability properties, but the entire shape of the valley 
region. 

For the S-types it is necessary that the cooling time- 
scale at every point of the cooling layer remains of the 
same order as the hydrodynamical time-scale. With the ex- 
ception of the rear boundary layer, the cooling layer must 
have temperatures associated with moderately positive (3. 
For the temperatures we encounter in the rear boundary 
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layer, the shape of the RLF may be arbitrary. In particu- 
lar, it can be strongly negative (e.g. for smooth types with 
post-shock temperatures in a valley region). This does not 
necessarily mean that the shape of the RLF for temper- 
atures in the rear boundary layer is unimportant for the 
properties of the instability. On the contrary, we notice the 
importance of the conditions in the rear boundary layer for 
the stability limit and the mode in many examples. 

For the modes the situation seems even more puzzling. 
We notice the absence of the overtone modes in the strong 
form of the F-modc oscillation. Wo ascribe this to the to- 
tal destruction of information associated with each col- 
lapse of the shell. Overtone modes just have no time to 
grow. We emphasize the existence of a strong form of the 
10-mode which is able to radiate a substantial amount of 
energy. Therefore, the 10-mode is potentially as impor- 
tant as the F-mode even though the temperature range 
of its appearance is considerably smaller. Finally, we find 
it remarkable that, despite the highly nonlinear character 
of radiative cooling, we find weak and strong forms of the 
overstability which, while appearing phenomenologically 
so differently, are governed by the same modes. 

Prom the insight into the working mechanism of the 
different types one may find constraints for the transitions 
between them. However, there is certainly no simple way 
to exactly determine the stability properties at every T^j . 

8.2. Temperature regions with different stability properties 

Based on the above discussion we propose five differ- 
ent stability regions for RLFs. In Fig. 17 we sketch a 
RLF which shows the typical behavior of most high- 
temperature RLFs (see e.g. Schmiitzlcr & Tscharmitcr 
1993 or GEC). We note again that the boundaries be- 
tween the different regions can be given only approxi- 
mately, since they cannot be derived by a simple rule. 

• Region I: There is always a low-temperature stable 
region, which is also present for fully time-dependent cool- 
ing. IGF and GEC found that shocks with speeds below 
about 130 km/s are stable. 

• Region II: Here we encounter smooth types of the 
overstability. In this region, we expect the shape of the 
entire cooling function to be important for the stability 
properties. In particular, the rear boundary layer attached 
to the cold, compressed sheet can play a significant role. 
The upper boundary of this region lies at a temperature 
which is slightly above the point where the logarithmic 
slope of the RLF becomes strongly negative. 

• Region III: We find strong runaway instabilities (C- 
or I-type, always F-mode) when T^^ lies in a region of 
the RLF with negative logarithmic slope and T^* is suf- 
ficiently above the low-temperature edge of the negative 
slope region. Under these conditions, the overstability is 
nearly independent of the shape of the RLF at lower tem- 
peratures. 
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Fig. 17. Sketch of a typical radiative loss function A(r) 
on a log-log scale. The solid line denotes the region we 
have extensively investigated, the dashed line denotes the 
high-tcmporaturc region we have not investigated. Shown are 
the five qualitatively different stability regions, region I: sta- 
ble; region II: smooth types; region III: runaway cooling; re- 
gion IV: valley region; region *: behavior not known; region V: 
Bremsstrahlung region. 

• Region IV: This region is associated with a valley in 

the RLF which determines the stability properties. Note, 
that the valley is larger than stability region IV (e.g. be- 
tween 1.6 • 10^ and 1 • 10^ K in RLF2). The shape of the 
RLF below the valley seems of minor importance. The 
strong form of the 10-mode (M-type) can only occur in 
such a region. In addition, we find in this region the strong 
form of the F-mode as well as the smooth form of the 10- 
mode and even stable behavior. For this region it seems 
hardly possible to predict whether or not, and in which 
mode and type, the overstability is present. Any real shock 
which is slightly disturbed can abruptly change into any 
of these types. In this sense, the shock reacts truly chaotic 
in this region. 

• Region V: High-temperature Bremsstrahlung region 
with /3 = 0.5. Although we have not done any experiment 
in this region we predict a smooth S2-type oscillation if we 
are sufficiently far away from the low-temperature bound- 
ary where the logarithmic slope flattens to (3 ^0.5. Near 
the low-tcmpcrature boundary the existence of a F-mode 
may be possible. Also in this region the rear boundary 
layer of the hot shell may be important for the stability 
property. 

8.3. Parameterized radiative loss function versus more 
elaborated models 

IGF and GEC emphasize the need for fully time- 
dependent radiation hydro dynamical calculations which 
explicitly include ionization, recombination and non- 
equilibrium line-cooling, and we completely agree with 
them. However, such computations are very costly. At 
present and even in the near future, systematic longtime 
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Radiative cooling overstability in colliding flows 


types 


characteristics 


F-mode 


lO-mode 


strength 


oscillation amplitude 


'''cool / 


collapsing 


cools temp, to 


secondary 








F-mode 


Tiydro 


shell 


nebular T 


shocks 


Sl-type 


S2-type 


weak 


order one 


« 1 always 


no 


no 


none 


I-type 




moderate 


one magnitude 


temp. <, 1 


slow collapse 


no 


weak 


C-type 


M-type 


strong 


up to two magnitudes 


temp. << 1 


yes 


yes 


strong 



Table 1. Schematic summary of the classiflcation and the characteristics for the spherically symmetric cooling overstability. 
T stands for temperature and temp, for temporarily. Note, the oscillation amplitude of the 10-mode is always of order one, also 
for the strong form. 



evolution and multi-dimensional computations cannot be 
done using such an approach. One has, therefore, to deal 
with the use of parameterized loss functions. Some of the 
major short comings of parameterized RLFs compared 
to fully time dependent computations have already been 
mentioned in Sect 2.2, and were discussed in some more 
detail by GEC and IGF. A key point is that the ionization 
structure of the plasma does not necessarily correspond to 
the gas temperature. The cooling history is often crucial 
for the determination of the ionization structure. Conse- 
quently, the correct cooling rate cannot be described by a 
function of the temperature alone as in the case of a RLE. 

In the runaway cooling range (ui„ ~ 200 km/s), GEC 
directly compare one example of the fully time-dependent 
model with simulations applying two differently parame- 
terized models (their Eig. 3). The global dynamics result- 
ing from the two different approaches qualitatively coin- 
cide. The oscillation period and amplitude differ by some 
percent. For the runaway cooling range this result is not 
surprising. However, in cooling regions II and IV we as- 
sume that the correct shock dynamics is less well mirrored 
when using a RLE. Here, the stability properties cannot 
be pinned down to the shape of the loss function around 
T^^. A time dependent computation in this regime may 
yield different overall dynamics. 

Generally, newly shocked matter is likely to be un- 
derionized, whereas gas on its cooling track tends to be 
overionized. In regions II and III underionization tends to 
stabilize the shock since the effective cooling function flat- 
tens compared to the parameterized function at this tem- 
perature (see Eig. 2 of GEC). The influence of a possible 
overionization at lower temperature will not change this 
picture, because lower temperatures are of minor impor- 
tance for runaway cooling. On the other hand, in the valley 
stability region IV, underionization tends to decrease the 
logarithmic slope of the effective loss function at 7^^, and 
has probably a destabilizing effect. Overionization along 
the cooling track may slightly reduce the size of the col- 
lapsing part at the rear end of the cooling layer. 

Partly based on the global similarity between fully 
time dependent computations and RLF computations 
GEC suggested that the overall stability properties of a 



radiative shock can be determined by the local (3 at the 
high temperature end of the recombination zone, at least 
if only small perturbations of the steady state shock are 
considered. As has been demonstrated in this paper, this 
local condition applies only in runaway cooling regions. 
Eor other stability types we assume that the importance 
of the cooling behavior at lower temperatures persists also 
in a time-dependent computation. 

Also GEC have shown that some important details are 
missed or are different in the parameterized model, even in 
cases where the overall dynamics arc relatively well cap- 
tured by parameterized cooling. In addition, ICE have 
found considerable differences between calculated spectra 
for equilibrium and non- equilibrium cooling shocks. In gen- 
eral, we feel that the question of how well the use of pa- 
rameterized radiative loss functions can replace fully time- 
dependent calculations badly needs further investigation. 
Not only as far as the global dynamics are concerned but 
also with regard to diagnostics. A future perspective may 
be the use of improved parameterized models, e.g. the use 
of several different radiative loss functions, to better ac- 
count for the different cooling histories 

9. Conclusion 

Before drawing conclusions from this systematic study of 
the long-term evolution of the radiative cooling instabil- 
ity in colliding flows we emphasize that these results hold 
for parameterized cooling of the form ■ A{T). Possible 
relevances for fully time-dependent cooling are discussed 
in Sect. 8. Even though most of the results are obtained 
by computing spherically symmetric flows, the conclusions 
are only reliable for planar flows, since the size of the ra- 
diative shock is much smaller than the radius of the struc- 
ture. 

We emphasize the many different manifestations of 
the cooling overstability. For both practically important 
modes, the fundamental and first overtone mode, we find a 
strong and a weak form of the overstability. The key point 
for the distinction between these two forms is the ratio of 
the two governing time scales, r] = "^cool/Tiydro' which 
in the strong forms may be temporarily very small but 
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which is always of order one in the weak forms. We, there- 
fore, propose the introduction of different phenomenolog- 
ical types as summarized in Table 1. 

The system always relaxes to the same asymptotic os- 
cillatory or stable solution which only depends on the 
stationary post-shock temperature. This asymptotic so- 
lution is independent of the sign and the size of the dis- 
turbance. It is possible to distinguish five qualitatively dif- 
ferent stability regions with different stability properties. 
Even though the importance of the different shapes of the 
RLF in different temperature regions can be roughly esti- 
mated, we find no easy criterion which determines the sta- 
bility properties, neither the stability limit, nor the mode, 
nor the type. In particular, wc find many examples in 
which the logarithmic slope of the RLF at the stationary 
post-shock temperature is not the only relevant parameter 
for the stability behavior. 

Probably not all radiative shocks are unstable above a 
certain threshold velocity of the stationary shock. Due to 
the existence of 'valleys' in the radiative loss function we 
find stable islands or smooth types of the instability even 
for high velocities. A necessary but not sufficient condition 
for their occurrence is a logarithmic slope of the radiative 
loss function which according to linear theory is stable 
at the stationary post-shock temperature. We also find a 
strong form of the first overtone mode, which exists in 
such valley regions. 

We find different evolutionary scenarios when radiative 
shocks run into a single density disturbance. The system 
reacts aperiodic only if the length and the mass of the 
disturbance are appropriately scaled. If they are too small 
we find only a moderate, short reaction but the system 
remains in the original limit cycle. If the length is too big 
the system has time to adjust itself and the oscillations 
are modulated according to the new pre-shock conditions. 
If the mass is too large, the shell collapses immediately 
and its history is destroyed. 

The cold compressed gas downstream of an interstel- 
lar radiative shock can be important for the stability of 
the shock. When smooth flows collide, this layer is to a 
minor degree (ten percent effects) dynamically important 
in the strong forms of the instability. However, a density 
disturbance in the flow on the order of a few percent of 
the mass of the layer leads to a significant dynamical influ- 
ence of the layer. As long as we encounter strong forms of 
the overstability, the cold dense layer shows a rich internal 
structure. 

The numerical method and the coarseness of the nu- 
merical grid can significantly influence the computed so- 
lutions. The entire instability or any particular mode can 
vanish only due to a too coarse mesh combined with an 
inappropriate numerical method. This effect is most prob- 
ably present in all multidimensional calculations. 
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Table 2. RLFl based on the loss function by Cook et al. 
(1989). 
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(4) 


4.06 
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6.9 


(10) 


2 
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4 
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0.50 
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Table 3. RLF2 based on the time-dependent loss function by 
Schmutzler and Tscharnuter (1993). 
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Model 


Wind bubble (WB) 




geometry 


characteristics 


wind of central star 


interstellar medium 




^min 


rmax 


A r 


RLF 


Tc/Th 


massloss 


vel. 


temp. 


N 


vel. 


temp. 




cm 


cm 


cm 




°kK 


Ms/t/ 


km/s 


°kK 


cm"^ 


km/s 


°kK 


WBl 


1.6 (17) 


1 (20) 


2.4 (13) 


1 


15/10 


1 (-4) 


1.6 (3) 


15 


1.5 





10 


WBle 


1.6 (17) 


1 (20) 


2.4 (13) 


1 


10/10 


1 (-4) 


1.6 (3) 


15 


1.5 





10 


WB2 


1.6 (17) 


1 (20) 


2.4 (13) 


2 


15/10 


1 (-4) 


1.6 (3) 


15 


1.5 





10 


WBlf 


1.6 (17) 


1 (19) 


3.0 (11) 


1 


15/10 


1 (-4) 


3.0 (3) 


15 


75 





10 


WB2f 


1.6 (17) 


1 (20) 


2.4 (13) 


2 


15/10 


1 (-4) 


2.8 (3) 


15 


15 





10 




Supernova remnant (SNR) 




geometry 


characteristics 


stellar blast 


interstellar medium 




'^min 


rmax 


A r 


RLF 


Tc/Th 


mass 


energy 


temp. 


N 


vel. 


temp. 




cm 


cm 


cm 




°kK 




ergs 


°kK 


cm"^ 


km/s 


°kK 


SNRl 


1 (16) 


1 (21) 


4.9 (14) 


1 


15/10 


50 


1 (51) 


15 


1 





10 


SNR2 


1(16) 


1(21) 


6.1 (13) 


2 


15/10 


50 


1 (51) 


15 


1 





10 




Colliding winds (CW) 




geometry 


characteristics 


present wind of CS 


earlier wind of CS 




"^min 


rmax 


A r 


RLF 




massloss 


vel. 


temp. 


massloss 


vel. 


temp. 




cm 


cm 


cm 




°kK 


Mq/j/ 


km/s 


°kK 


Mq/j/ 


km/s 


°kK 


CW 


5 (15) 


1 (17) 


2.3 (10) 


2 


15/15 


1.1 (-7) 


2.5 (3) 


15 


2 (-6) 


20 


15 




Colliding planar flows (CPF) 




Instability 


geometry 


characteristics 


driving wind 


interstellar medium 




type 


domain 


A r 


RLF 


Tc/Th 


density 


vel. 


temp. 


density 


vel. 


temp. 






cm 


cm 




°kK 


cm"^ 


km / s 


°kK 


cm'"^ 


km/s 


°kK 


CPF_25 


C 


1 (20) 


2.4(13) 


2 


15/10 


0.025 


2800 


15 


15 





10 


CPF_35 


M 


1 (20) 


2.4(13) 


2 


15/10 


0.035 


2800 


15 


15 





10 


CPF.45 


S2 


1 (20) 


2.4(13) 


2 


15/10 


0.045 


2800 


15 


15 





10 


CPF_55 


stable 


1 (20) 


2.4(13) 


2 


15/10 


0.055 


2800 


15 


15 





10 



Table 4. Sample of computed colliding smooth flows. As a representative for the wind bubble flow we take NGC 6880 (WR 136). 
We assume model parameters as derived by Schmutz et al. (1989). In addition, we compute models with faster winds (models .f). 
For the supernova remnant we assume the same parameters as Bertschinger (1986). 80 percent of the initial energy is thermal in 
our models. For the colliding wind we take the same parameters as Icke et al. (1992) took for their (2D) model of an (aspherical) 
planetary nebulae for the colliding wind model, tj^^^ and rmax denote the minimum and maximum extent of the computational 
domain. Ar denotes the finest spatial discretization of the computational mesh. Tq denotes the cooling limit, Tjj the heating 
limit. CS stands for central star. 
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Table 5. Sample of the disturbed flows discussed in Sect. 6. We disturb model WBl (see Table 4) according to Eq. 13 at a time 
of approximately 69'000 years where we encounter typical I- type overst ability. The period of the limit cycle, tq, is approximately 
1100 years. The mean mass Mjg of the leading shell is approximately 6 • 10^^ gram. The rear boundary of the leading shell has a 
mean velocity of 9.95 ■ lO'' cm/s. Multiplication with the period of the limit cycle yields a characteristic length Ic ~ 3.4- lO^'^ cm. 
The mass M^^j of the cold dense layer is 7.1 • 10^® gram. denotes the mass, Ao the length and § the amplitude of the 
disturbance. 0o denotes the phase of the limit cycle at which the leading shock hits the disturbance. Trei denotes the relaxation 
time for the aperiodic cases {ap). The other cases are: nop: noisy original periodic evolution, mp: modulated periodic evolution, 
(co;) collapse of the leading shell. Note that some models may be listed several times to better illustrate dependences. 



Plewa T., 1995, MNRAS 275, 143 
Ryu D., Vishniac E. T., 1991, ApJ 368, 411 
Schmutz W., Hamann W.-R., Wessolowski U., 1989, A&A 
210, 236 

Schmutzler T., Tscharnuter W. M., 1993, A&A 273, 318 
Strickland R., Blondin J. M., 1995, ApJ 449, 727 
Waldcr R., 1993, Ph.D. thesis, ETH Ziirich 
Wolff M. T., Gardner J. H., Wood K. S., 1989, ApJ 346, 
833 



This article was processed by the author using Springer- Verlag 
A&A style file L-AA version 3. 



